Chapter 5 



Methods for ordinary 
differential equations 



5.1 Initial- value problems 



Initial-value problems (IVP) are those for which the solution is entirely known 
at some time, say t — 0, and the question is to solve the ODE 



for other times, say t > 0. We will consider a scalar y, but considering 
systems of ODE is a straightforward extension for what we do in this chapter. 
We'll treat both theoretical questions of existence and uniqueness, as well as 
practical questions concerning numerical solvers. 

We speak of t as being time, because that's usually the physical context 
in which IVP arise, but it could very well be a space variable. 

Does a solution exist, is it unique, and does it tend to infinity (blow up) 
in finite time? These questions are not merely pedantic. As we now show 
with two examples, things can go wrong very quickly if we posit the wrong 



y'(t) = f(t,y(t)), 



y(o) = yo, 



ODE. 



Example 12. Consider 

V = Vy, 

By separation of variables, we find 



y(0) = 0. 





it + Cf 
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Imposing the initial condition yields C = 0, hence y(t) = t 2 /4. However, 
y(t) = is clearly another solution, so we have non-uniqueness. In fact, there 
is an infinite number of solutions, corresponding to y(t) = for < t < t* for 
some t* , which then takes off along the parabola the parabola y(t) = (t — t*) 2 /4 
for times t >t*. 

Example 13. Consider 

y' = y 2 , y(0) = 1. 
By separation of variables, we find 

The initial condition gives C = —1, hence y(t) = 1/(1 — t). It blows up at 
time t = 1, because y(t) has a vertical asymptote. We say that the solution 
exists locally in any interval to the left of t = 1, but we don't have global 
existence. 

Blowup and non-uniqueness are generally, although not always 1 , unreal- 
istic in applications. A theory of existence and uniqueness, including global 
existence (non blowup), is desired to guide us in formulating valid models 
for physical phenomena. 

The basic result is the following. 

Theorem 8. (Picard) For given T,C, and y , consider the box B in (t,y) 
space, given by 

B = [0,T\ x [y Q -C,y Q + C\. 

Assume that 

• f(t,y) is continuous over B; 

• \f(t,y)\ < K when (t,y) G B; (boundedness) 

• \f(t,u) — f{t,v)\ < L\u — v\ when (t,u),(t,v) G B. (Lipschitz 
continuity). 

1 ln nonlinear optics for instance, laser pulses may "collapse". This situation is some- 
what realistically modeled by an ODE that blows up, although not for times arbitrarily 
close to the singularity. 
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Assume furthermore that C > j;(e LT — 1). Then there exists a unique y E 
C^T], such that 

y'(t) = f(t,y(t)), 2/(0) =2/0, 

and such that \y(t) — y \ < C . In short, the solution exists, is unique, and 
stays in the box for times < t < T . 

Proof. The technique is called Picard's iteration. See p. 311 in Suli-Mayers. 

□ 

The ODE y' = ^fy does not satisfy the assumptions of the theorem above 
because the square root is not a Lipschitz function. The ODE y' = y 2 does 
not satisfy the assumptions of the theorem above because y 2 is not bounded 
by a constant K. 

5.2 Numerical methods for Initial- Value Prob- 
lems 

Here is an overview of some of the most popular numerical methods for 
solving ODEs. Let t n = nh, and denote by y n the approximation of y(t n ). 

1. Forward Euler (a.k.a. explicit Euler). 

Vn+i = Vn + hf(t n ,y n ). 

This formula comes from approximating the derivative y' at t — t n by 
a forward difference. It allows to march in time from the knowledge of 
y n , to get y n+1 . 

2. Backward Euler (a.k.a. implicit Euler). 

y n+ i = y n + hf(t n+1 ,y n+1 ). 

This time we use a backward difference for approximating the derivative 
at t — t n+ i. The unknown y n+1 appears implicitly in this equation, 
hence the name implicit. It still needs to be solved for as a function 
of y n , using (for instance) Newton's method. The strategy is still to 
march in time, but at every step there is a nonlinear equation to solve. 
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3. Trapezoidal (a.k.a. midpoint) rule (implicit). 

y n +i =y n + ^ [f(t n , y n ) + /(t„+i, y n +i)} ■ 

In this equation Vn+1 ~ yn has the interpretation of a centered difference 
about the midpoint t n+ i = fn+ *" +1 , but since f(t n+ i,y n+ i) is not ac- 
cessible (y n+ i is not part of what we wish to solve for), we replace it 
by the average \ [f(t n , y n ) + f(t n+ i, y n +i)]- This gives a more balanced 
estimate of the slope Vn+1 ~ Vn , It is an implicit method: y n+ i needs to 
be solved for. 

4. Improved Euler, Runge-Kutta 2 (explicit). 

y n +i = y n + hf(t n ,y n ), 

y n +i =y n + ^ [f(t n , y n ) + /(* n+ i, y n +i)] ■ 

This is the simplest of "predictor-corrector" methods. It is like the 
midpoint rule, except that we use a guess y n +\ for the unknown value 
of y n +i in the right-hand side, and this guess comes from the explicit 
Euler method. Now y n+ \ only appears in the left-hand side, so this is 
an explicit method. 

5. Runge-Kutta 4 (explicit). 

y n +i = y n + h[h + 2k 2 + 2k 3 + fc 4 ], 
where the slopes k±, . . . , k± are given in succession by 

h = f(tn,y n ), h = f(t n + ^,y n + ^h), 

h = f(t n +^,y n + ^h), h = f(t n + h,y n + hk 3 ). 

6. There are also methods that involve not just the past value y n , but a 
larger chunk of history y n _i, y n _ 2 ,etc. These methods are called multi- 
step. They are in general less flexible than the one-step methods de- 
scribed so far, in that they require a constant step h as we march in 
time. 
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Two features of a numerical method are important when choosing a nu- 
merical method: 

• Is it convergent, i.e., does the computed solution tend to the true solu- 
tion as h — > 0, and at which rate? 

• Is it stable, i.e., if we solve with different initial conditions y and y , 
are the computed solutions close in the sense that \y n — y n \ < C\y — y \, 
with n = 0(1 /h), and C independent of hi 

5.2.1 Convergence 

To understand convergence better in the case of one-step methods, let us 
write the numerical scheme as 

y n+1 = V(t n ,y n ,h), 

and introduce the local error 

e n+1 (h) = ty(t n ,y(t n ),h) -y(t n+1 ), 

as well as the global error 

E n {h) =y n -y(t n )- 

The expression of the local error can be explained as "trying the exact 
solution in the numerical scheme" — although the exact solution is unknown. 
It is a lot easier to approach the convergence question via local errors than 
global errors. It would not be practical, for instance, to "try an approximate 
solution in the exact ODE". 

Convergence is the study of the global error. Consistency is the study 
of the local error. A numerical method is called consistent if the local error 
decays sufficiently fast as h — > that there is hope that the global error 
would be small as well. The particular rate at which the local error decays 
is related to the notion of order of an ODE solver. 

Definition 6. (Consistency) ^ is consistent if, for any n>0, 

lim^ = 

h^O h 
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Definition 7. (Order) \P is of order p if e n (h) = 0(h p+1 ). 

The basic convergence theorem for one-step solvers, that we 
will not prove, is that if the local error is 0(h p+1 ), then the global 
error is 0(h p ). This convergence result is only true as stated for one- 
step methods; for multi-step methods we would also need an assumption of 
stability (discussed below). Intuitively, the local errors compound over the 
0(1/ h) time steps necessary to reach a given fixed time t, hence the loss of 
one power of h. Of course the local errors don't exactly add up; but they do 
up to a multiplicative constant. It is the behavior of the global error that 
dictates the notion of order of the numerical scheme. 

It is a good exercise to show, using elementary Taylor expansions, that the 
explicit and implicit Euler methods are of order 1, and that the midpoint rule 
and improved Euler methods are of order 2. It turns out that Runge-Kutta 
4 is of order 4, but it is not much fun to prove that. 



5.2.2 Stability 

Consistency and convergence do not tell the whole story. They are helpful 
in the limit h — > 0, but don't say much about the behavior of a solver in the 
interesting regime when h is small, but not so small as to make the method 
computationally inefficient. 

The single most important consideration in the regime of moderately 
small h is perhaps stability. The criterion stated above for stability (\y n — 
y n \ < C\yo — yo\) is too complex to handle as such. The important ideas 
already appear if we study the representative setting of linear stability. The 
linearization of the ODE y' = f(t, y) about a point yo is obtained from writing 
the Taylor expansion 

j t (y(t)-y ) = f(t,y(t)) = f(t,y ) + ^-(t, y )(y(t) - y ) + o(\y(t) - y \). 

The culprit for explosive (exponential) growth or decay is the linear term 
T^(t,y )(y(t) — y ). (Indeed, if the other two terms are neglected, we can 
write the solution, locally, as an exponential.) For practical purposes it is 
sufficient to check stability for the linear equation y' = Xy, keeping in mind 
that A is a number representative of the derivative §^(i, yo) of / at yo in the 
y variable. 
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Definition 8. (Linear stability 2 ) Suppose y' = Xy for some A G C. Then the 
numerical method ^ is linearly stable if y n — ► as n — > oo. 

Of course linear stability depends on the value of A. Stability for the 
original equation y' = Xy is guaranteed if Re(A) < (because the solution is 
y(0)e xt ), and the question is that of showing whether a numerical method ^ 
is stable under the same condition or not. 

If a numerical method is stable in the above sense for a certain range of 
values of A, then it is possible to show that it will be stable for the ODE 
y' — f(t, y) as long as ^ is in that range of A (and / is smooth enough). We 
won't prove this theorem here. 

Let us consider a few examples 

Example 14. For the forward Euler method applied to y' = Xy, we get 

y n+ i = y n + hXy n = (1 + hX)y n . 

The iterates y n tend to zero provided |1 + hX\ < 1, where the \ ■ \ denote 
the complex modulus. Write hX — x + iy, so that the inequality becomes 
(1 + x) 2 + y 2 < 1. This is the equation of a disc in the complex (hX)-plane, 
with center at — 1, and radius 1. If hX sits inside this disk, the method 
is stable, and otherwise it isn't. We say that the forward Euler method is 
conditionally stable: typically, we require both Re(X) < and a small step 
size h in order to guarantee stability. 

Example 15. The backward Euler method applied to y' = Xy gives 

y n +i = y n + hXy n+ i, 

or in other words 

y n 

The iterates y n tend to zero provided |1 — hX\ > 1. In terms of hX = x + iy, 
this becomes (x — l) 2 + y 2 > 1. This condition is satisfied whenever hX is 
outside the disc in the complex (hX)-plane with center at +1, and radius 1. 
In that case the method is stable, and otherwise it isn't. We say that the 
backward Euler method is unconditionally stable: the stability zone for the 
ODE (Re(X) < 0) is always included in the stability zone of the numerical 

2 Sometimes called A-stability in some texts. 
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< 1. 



method, regardless of h. In fact the zone of stability for the backward Euler 
method is larger than the left complex half-plane, so there exist choices of 
(large) time steps for which the numerical method is stable although the ODE 
isn't. 

Example 16. The linearized stability analysis of the midpoint method gives 

Vn+l =Vn + h\{— + —^-h 

hence 

f l + hX/2 \ 
Vn+i- [ 1 _ hX / 2 ) yn - 

The method is stable provided 

1 + hX/2 
1 1 - h\/2 1 

In terms of hX = x + iy, we can simplify to get 

(1 + |) 2 + (|) 2 <(1-|) 2 + (|) 2 , 

which is true if and only if x < 0. So the stability region is Re(hX) < 0, the 
same as that of the ODE. As a result, the method is unconditionally stable. 

It is a general rule that explicit methods have conditional stability (sta- 
bility only happens when the time step is small enough, if it does at all), 
whereas implicit methods are unconditionally stable. 

The stability regions for the Runge-Kutta methods are plotted on page 
351 of Suli-Mayers. 

Stability can also be studied for systems of ODEs y'(t) = f(t, y(t)) where 
both y and / are vectors. The interesting object is now the Jacobian matrix 

A = V y f(t,y ). 

(it is a matrix because the gradient of a vector function is a "vector of vec- 
tors", i.e., a matrix.) The linearized problem is now 

y'(t) = Ay(t), 

which has the solution y(0)e tA , where e tA is a matrix exponential. Stability 
will happen when the eigenvalues A of A all obey Re(A) < 0. So the story 
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is the same as in the scalar case, except that A is now any eigenvalue of 
the Jacobian matrix. So we need to make sure that hX is in the stability 
zone of the ODE solver in the complex plane, for each A an eigenvalue of the 
Jacobian matrix. 

Problems for which A has a very large, negative real part are called stiff. 
Physically they are very stable, but they pose numerical problems for ex- 
plicit methods since the region of stability does not extend very far along 
the negative real axis. Implicit methods are hands down the best for stiff 
problems. 

5.2.3 Miscellaneous 

Another interesting class of methods, and a way to naturally design high- 
order methods, is deferred correction. Assume that time is split into a uniform 
grid tj = jh, and that some low-order method has given us the samples 
yi, ... , y n+ i for some (small) m. Denote by n n (t) the n-th order interpolation 
polynomial passing through the (tj,yj). Consider the error ("defect") 

S(t)=y(t)-Tr n (t). 

It obeys the equation 

5'(t) = f(tMt))-<(t), m = o. 

We do not have access to y(t) in the argument of /, but we can replace it by 
our best guess n n (t) +S(t). This way we can compute an approximate defect 

~5'(t) = f(t,7r n (t) + ~5(t))-7r' n (t), 6(0) = 0. 

using the same (or another) low-order method. Once S(t) is computed, add 
it back to the interpolant to get 

y(t) = n n {t) + ~5{t). 

This procedure can be repeated a few times to get the desired accuracy. 
(This version of deferred correction is actually quite recent and due to Dutt, 
Greengard, and Rokhlin, 2000). 
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5.3 Boundary- value problems 

Boundary-value problems (BVP) are ODE where some feature of the solution 
is specified at two ends of an interval. The number of initial or boundary 
conditions matches the order of the highest derivative in the equation, hence 
such ODE are generally second-order scalar equations. The simplest exam- 
ples are 

-u"{x) = f(x), x e [0, 1], u(0) = a,u(l) = b (Dirichlet) 

-u"(x) = f(x), x G [0, 1], u'(0) = a,u'(l) = b (Neumann) 

-u"(x) = f(x), x e [0, 1], u(0) = u(l) (periodic) 

We don't really mean to study the equation —u" = f for its own sake (you can 
find the general solution by integrating / twice and fixing the two constants 
from the boundary conditions), but its study serves two purposes: 

• It is the simplest ODE that models problems in continuous elasticity: 
here u is the displacement of a vertical elastic bar, or rod, that sags 
under its own weight. The right-hand-side / is the gravitational force 
as a function of x, and u' is the "elongation", or strain of the bar. 
A condition where u is specified means that the bar is fixed at that 
end, while the condition u' — would mean that that end is free. The 
condition u(0) = u(l) means that it's an elastic band. By solving the 
ODE we are finding the displacement u generated by the force /. 

• It is the simplest boundary-value problem to treat numerically, and 
contains many of the important features of such problems. It needs to 
be understood before moving on to any other example. 

Alternative problems of this kind are 

—u"(x) + a(x)u(x) = f(x), u(0) = a, u(l) = b, 

for instance, the solution of which does not generally obey an explicit formula. 

A first intuitive method for solving BVP is the shooting method. Consider 
again u"(x) + a(x)u(x) = f(x), u(0) = a, u(l) = b. We cannot in general 
march in time from to 1 (the variable x is more often a spatial variable), 
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but we can guess what the derivative should have been for the solution to 
end up near b at x — 1. Introduce a parameter s, and write 

— u"(x; s) + a(x)u(x; s) = f(x), u(0] t) = a, u'(0] s) = s. 

Of course in general we won't reach b at x = 1, but we can refine our estimate 
of s so that we get closer to it. The problem becomes that of solving for s in 
u(l; s) = b, where u(l; s) is defined implicitly from the ODE. This equation 
that can be viewed as a root-finding or fixed-point problem and solved using 
any of the methods we've seen previously. The secant method only requires 
evaluations of u(l; s), which involves solving the ODE. To set up Newton's 
method, we need the derivative of the solution as a function of its 

parameter s. Exercise: find the value of this derivative by differentiating the 
ODE in the t parameter, and obtain an ODE for ^(1; s) itself. 

The shooting method is easy and intuitive (and there isn't much more 
to say about it), but it has the disadvantage of being restricted to ODE. In 
contrast, we'll now investigate finite difference methods, which can also be 
applied to partial differential equations (PDE). 

Let's start by considering the problem — u" = f, u(0) = u(l) = 0. Con- 
sider the grid Xj = jh, j — 0, . . . , N, h — l/N. This grid has N + 1 points. 
The usual 3-point stencil for the centered second difference gives rise to 



U. 



j+i — 2Uj + Uj-i , . 
" u l — = f( x j)- 



(5.1) 



In this context, capital U denotes the numerical solution, while lowercase 
u denotes the exact solution. For the boundary conditions, we simply set 

. N — 1 are unknowns and 
F generated by (5.1). The 



Uq = Un = 0. The rest of the Uj for j — 1, . 
can be solved for from the linear system KU = 
resulting matrix K (of size iV — 1 by iV — 1) is 



(2 -1 
-1 2 -1 



K = 



h 2 



\ 



-1 2 -1 
-1 2 



The zero elements are not shown. The right-hand side is here Fj = f(xj). 
In Matlab, one then solves U as K\F. Had the boundary conditions been 
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Uq = a and Un = b instead, it is a good exercise to check that this does not 
change K, but that the right-hand side gets modified as 



F = 



\ 



f(x N -2) 
\f(x N -l) + J?) 

Of course, the matrix K should be invertible for this strategy to make 
sense. We will see below that this is the case. It is also important that K~ l 
be bounded for convergence of the numerical scheme. 

Definition 9. The local truncation error (LTE) of a numerical scheme KU = 
F, is the error made when evaluating the numerical scheme with the exact 
solution u(xj) in place of the numerical solution Uj. It is the quantity r in 



Ku = F + r. 



The local truncation is directly obtained from the truncation error of the 
finite-difference scheme, for instance 

_ ^ + i)-2^)+ M (x 3 -_ 1 ) = fM + Q{h% 
so the LTE is 0(h 2 ). 

Definition 10. The (actual) error of a numerical scheme KU = F , is the 
vector of differences ej = u(xj) — Uj. 

In order to obtain the error e from the LTE, one writes 

Ku = F + t, KU = F, 

and subtract those two equations. This gives K(u — U) = (F — F) + r, or 
Ke = t. If K is invertible, this gives 

e = K~ x t. 



The next few sections introduce the tools needed to control how large e 
can get from the knowledge that it is r "magnified" by K v . 
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5.3.1 Matrix norms and eigenvalues 

This section is mostly a linear algebra refresher. 

Definition 11. The spectral norm, or 2-norm, of (any rectangular, real) 
matrix A is 

\\M\2 

LA 2 = max ",, „" , 

\\x\\2 

where the vector 2-norm is \\x\\ 2 = x 1- ^ n other words, the matrix 2- 

norm is the maximum stretch factor for the length of a vector after applying 
the matrix to it. 

Remark that, by definition, 

||Ae|| 2 < PH2IMI2, 

for any vector x, so the matrix 2-norm is a very useful tool to write all sorts 
of inequalities involving matrices. 

We can now characterize the 2-norm of a symmetric matrix as a function 
of its eigenvalues. The eigenvalue decomposition of a matrix A is, in matrix 
form, the equation AV = VA, where V contains the eigenvectors as columns, 
and A contains the eigenvalues of the diagonal (and zero elsewhere.) For those 
(non-defective) matrices for which there is a full count of eigenvectors, we 
also have 

A = VAV' 1 . 

Symmetric matrices have a full set of orthogonal eigenvectors, so we can 
further write V T V = J, VV T = I (i.e. V is unitary, a.k.a. a rotation), hence 

A = VAV T . 

In terms of vectors, this becomes 

A = ^ VjXivf. 

i 

Note that \ are automatically real when A = A T . 
Theorem 9. Let A = A T . Then 

P|| 2 = max|A,(A)|. 
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Proof. First observe that the vector 2-norm is invariant under rotations: if 
V T V = I and VV T = I, then ||Va;||2 = 1Mb- This follows from the defini- 
tion: 

\\Vx\\l = (Vx) T Vx = x T V T Vx = x T x = \\x\\l 
Now fix x and consider now the ratio 
\\Ax\\ 2 2 \\VAV T x\\ 2 2 



I-*' || 2 IKII2 

\\AV T x\ 



1 1 *^ 1 1 2 

IIAyjll 
\\Vy\\l 
JJAyJH 

llvlli 

2 

2 Hi 



(unitary invariance) 
(change variables) 
(unitary invariance again) 



E*y 



This quantity is maximized when y is concentrated to have nonzero compo- 
nent where Aj is the largest (in absolute value): yj = 1 when |Aj| = max n |A n |, 
and zero otherwise. In that case, 

' I,, = max A„ , 
Iklll 

the desired conclusion. □ 

Note: if the matrix is not symmetric, its 2-norm is on general not its 
largest eigenvalue (in absolute value). Instead, the 2-norm is the largest 
singular value (not material for this class, but very important concept.) 

One very useful property of eigenvalues is that if A = VAV T is invertible, 
then 

A' 1 = VA' l V T 
(which can be checked directly), and more generally 

f(A) = Vf{K)V T } 

where the function / is applied componentwise to the \. The eigenvectors of 
the function of a matrix are unchanged, and the eigenvalues are the function 
of the original eigenvalues. 
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If a matrix A is not invertible, things usually go wrong when trying to 
solve the linear system Ax = b. 

Definition 12. The nullspace of (any rectangular, real) matrix A is the space 
Null(A) of all vectors v such that Av = 0. In other words, it is the eigenspace 
corresponding to the eigenvalue zero. 

Null (A) always contains the zero vector. The following conditions are 
equivalent to characterize singular matrices: 

• A is singular (non-invertible); 

• Null (A) contains some nonzero vector; 

• is an eigenvalue of A; 

• det(A) = 0; 

• The rows/columns are linearly dependent; 

• (Zero is a pivot in the row echelon reduction of A.) 

We now present a version of the inversion theorem for symmetric matrices. 
If the matrix is not symmetric, the statement looks quite different. 

Theorem 10. Let A = A T . Consider the system Ax = b. 

• If A is invertible, then the solution is unique and x = A~ 1 b. 

• If Null(A) = span(vi, . . . v m ) ^ ; then 

— If b has a component along any of the Vj (i.e., vjb ^ for some 
j), then the system has no solution. 

— If all vjb — 0, j — 1, . . . , m, then there exists an infinite number of 
solution to the system. Ifx is a solution, then so is ^o + X^jli c j v j 
for arbitrary coefficients Cj. 

In terms of eigenvectors Vj, if the matrix is invertible, the solution of 
Ax = b is 



j'=i 
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If A is not invertible, but vjb = for all the eigenvectors Vj, j — 1, . . . , m 
corresponding to the zero eigenvalue (as in the theorem above), then we still 
have 

N m 
j=m+l ^ j=l 

where the first sum only runs from m + 1 to N, and the coefficients Cj are 
arbitrary. (Apply the matrix A to this equation to see that they disappear.) 

If vjb 7^ 0, then the operation j^vjb would result in an infinity when 
Xj = — a crude reason to see why the system has no solution in that case. 



5.3.2 Properties of the matrix K 

We are now equipped to study the matrix K and its inverse. Recall that we 
are interested in controlling e = K~ x t to get the error from the LTE. From 
the results in the previous section, we know that 

||e|| 2 = \\K- x t\\ 2 < WK-^Urh. 



Since K is a symmetric matrix, 



1 



\\K'% = maxlA^TT- 1 )! = . — — ■ 
J mm.lA^fOI 

So it remains for us to understand the eigenvalues of K, and specifically to 
show that the minimum eigenvalue (in absolute value) does not get too small. 
What we mean by this is that there should exist a number c > such that 

c < min \ Xj(K)\, 
j 

independently of the grid spacing h (recall that the size of K and its entries 
depend on h.) Only in that case will we be able to conclude that ||e|| is of 
the same order of magnitude as ||r||. 

Note in passing that if r = 0(h 2 ) then ||r|| 2 = \/Ei T i wil1 be °(\f\^) = 
0(h 3 / 2 ), which is not very appealing. Instead, it is common to modify the 
vector 2-norm as 



Mb, - ^E* 
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so as to restore ||t||2,/i = 0(h 2 ) (note how h times the sum resembles an 
integral quadrature.) In that case, we also expect ||t||2,/i = 0(h 2 ). The 
reasoning with matrix 2-norms does not change one bit form this different 
choice of normalization. 

Let us now study the eigenvalues and eigenvectors of K. The best way 
to guess them is to notice that KU = F is a discretization of — u" = f with 
Dirichlet boundary conditions. The eigenvalue problem 

-v" = Xv, v(0) = v(l) = 0, 

has a solution in terms of sines: for each n > 1, we have the pair 

v n (%) — sin(n7rx), A = n 2 n 2 . 

This analogy is very fruitful: the eigenvectors of K are precisely the v n 
sampled at Xj = jh, 

vf ] = sm(nirjh), j = 1,...,N -1, n = l,...,N -1. 

(here n is the label index, and j is the component index.) It is straightforward 
and a little tedious to check from trigonometric formulas that v ^ defined by 
this formula are indeed eigenvectors, with eigenvalues 

4 . 2/ Tmh. 
A n = — sm 2 (— ), n = l,...,N-l. 

A Taylor expansion for small h shows that the minimum eigenvalue is Ai = 
7r 2 + 0(h 2 ), and this 0(h 2 ) is positive, so that Ai > n 2 . 

Note in passing that since K is symmetric, the eigenvectors are automat- 
ically orthogonal, hence there is a way to normalize them so that V T V = I 
and VV T = I. Applying the matrix V T is called the discrete sine trans- 
form (DST), and applying the matrix V is called the inverse discrete sine 
transform (IDST). 

The formula for the eigenvalues has two important consequences: 

• The matrix K is invertible, now that it is clear that all its eigenvalues 
are positive. So setting up the discrete problem as KU = F makes 
sense. 

• Returning to the error bound, we can now conclude that \\e\\2,h < 
\ IMUh = 0(h 2 ), establishing once and for all that the finite-difference 
method is second-order accurate for — u" = f with Dirichlet boundary 
conditions. 
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The reasoning of first establishing consistency (small LTE) and showing 
a stability result transferring at the level of the actual error is very common 
in numerical analysis. Consistency + Stability = Convergence. 

5.3.3 Other boundary conditions and other equations 

The boundary conditions (BC) play a crucial role in the description of the 
problem: if the BC changes, the solution changes completely, and so can the 
LTE and the error. For instance, let's return to the Neumann problem 

-u"(x) = f(x), x G [0, 1], «'(0) = a,u'(l) = b (Neumann) 

The discretization of the interior equations is the same as previously, but at 
least one extra equation needs to be added for the BC. For instance, at zero, 
we may write a forward difference 

C/i-C/o 

and similarly, a backward difference at x — 1. Each BC adds one row and 
one column to the original matrix K, and results in a different system of 
equations. The choice above leads to a first-order LTE, only 0(h), even 
though the LTE for the interior equations is 0(h 2 ). This is enough to spoil 
the error itself: we don't have ||e|| 2 ,h = 0(h 2 ) in general, as a result of the 
low-accuracy boundary conditions. 

A more accurate way to discretize a Neumann boundary condition is to 
introduce a "ghost" node U-i, and write 

Ui - U-i 
-^h— = a - 

This new unknown is linked to the others by writing the equation one more 
time at j — 0, 

= U 1 + 2U -U„ l 

p = fM- 

(Previously, we had only evaluated the equation at j — 1, . . . , ./V — 1.) This 
results in two additional rows and columns per BC. The same treatment 
should be applied at x — 1. 

Note that these additional "boundary" rows can be scaled so that the 
resulting matrix looks symmetric, or at least more balanced. For instance, 
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a by 2/h (including the right-hand side) we get 

f-1 1 \ 

-1 2 -1 

-1 2 -1 

v 

The eigenvalues depend on this choice of normalization! 

Notice that, unlike the Dirichlet problem, the Neumann problem has a 
nullspace. The vector identically equal to 1 is in the nullspace of either of the 
2 discretizations we have presented. As a result, there exists a solution only 
if the admissibility condition l T f = ^ifixi) = is satisfied (see a theorem 
in the previous section). This is also a feature of the non-discretized BVP: 
it has a solution if and only if f(x) dx = 0. 

The periodic problem is also very interesting: 

-u"(x) = f(x), x G [0, 1], u(0) = u(l) (periodic) 

The boundary condition U — Un modifies K by adding two elements in the 
bottom-left and top-right: 

(2 -1 -1\ 

-1 2 -1 

-1 2 -1 
\-l -1 2 J 

This new matrix C is circulant and singular. We can check that the vector 1 is 
in its nullspace, so the BVP has a solution if and only if l T f = £V f{xi) = 0. 
The eigenvectors of C are the Fourier components 

Since C is symmetric, these v <ra * ) are orthogonal (and orthonormal when di- 
vided by y/~N). To deal with vectors and matrices that have complex entries, 
don't forget that the transposes come with a complex conjugate, so the dot 
product is 

x*y = x T y = ^XiVi. 

i 

The norm is ||x||2 = \ x %\ 2 - The orthogonality relations for the matrix of 

eigenvectors are now V*V = I and VV* = I, where * is transpose conjugate. 



by rescaling the 
the resulting matrix 



T 
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